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SUMMARY 


A description of a digital computer simulation of an Annular Momentum 
Control Device (AMCD) laboratory model is presented. The AMCD is a momentum 
exchange device which is under development as an advanced control effector for 
spacecraft attitude control systems. 

The digital computer simulation of this device incorporates the following 
models: six degree-of-freedom rigid body dynamics, rim warp, controller 

dynamics, nonlinear distributed element axial bearings , as well as power driver 
and power supply current limits. An annotated FORTRAN IV source code listing 
of the computer program is included. 


INTRODUCTION 

Momentum exchange devices are in use on a number of spacecraft and have 
demonstrated their ability to perform attitude control reliably and efficiently 
for vehicles as large as SKYLAB. The Annular Momentum Control Device (AMCD) 
concept was an outgrowth of research at the NASA Langley Research Center 
(NASA-LARC) directed toward development of momentum exchange devices. Reference 
1 presents a detailed description of the AMCD concept as well as potential 
applications. 

In order to investigate any potential problems in implementing the AMCD 
concept, a laboratory model of the AMCD was designed and fabricated under 
contract. The laboratory model has been delivered to NASA-LARC and preliminary 
tests of the device have been performed. A detailed description of the 
laboratory model is presented in reference 2. Reference 3 presents the results 
of the static and low-speed dynamic tests which include spin motor torque 
characteristics as well as spin motor and magnetic bearing drag losses. 

A digital computer simulation of the laboratory model AMCD has been 
developed as an analytical tool to investigate implementation problems with the 
laboratory model. This report presents a description of that simulation. The 
simulation incorporates six degree-of-freedom rigid body rim dynamics with rim 
warp superimposed, as well as power driver and power supply current limits, 
Hquations of motion of the rim are not developed in this report but are taken 
from reference 4. Similarly, the mathematical model of the nonlinear axial 
bearing element is taken from reference 5. 
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This description includes source code for the simulation computer program 
which is written in FORTRAN IV. An annotated listing of the program is 
contained in appendix A. 


SYMBOLS 


CNG 

CPS 

CR 

CRBAR 

CRL 

CRLIM 

CRSM 

CRU 

DEG 

DELT 

DELTPR 

D1V1 , DIV2 

FA, FB, FC 

FA1 , FA2 , FA3, FA4 

FBI, FB2, FB3 , FB4 

FC1 , FC2 , FC3, FC4 

FCNY 

FRA, FRB, FRC 

FI, F2 , F3 

GAA, GAB, GAC 

GAAC, GABC, GACC 

GAACP , GABCP, GACCP 

GAA1 , GAA2 , GAA3 , GAA4 

GAB1 , GAB2, GAB3 , GAB4 

GAC1 , GAC 2 , GAC3, GAC4 

GPBAR 

GPL 

GPPP 

GPRMS 

GPSM 

GPSQSM 

GPU 

GRA, GRB, GRC 
GRAP , GRBP, GRCP 
HI, H2, H3 
H1DT, H2DT, H3DT 
HI DTP, H2DTP, H3DTP 
IAA, IAB, I AC 
IBIAS 

INA1 , INA2 
IN A3 

IRA, IRB, IRC 
KAA, KAR 
KRA, KRR 


negative current limit hit counter 
positive current limit hit counter 
instantaneous power supply current 
average power supply current 
interval minimum power supply current 
power supply current limit 
interval sum of power supply current 
interval maximum power supply current 
radian to degree conversion factor 
integration step size 
output print interval 

temporary variables in data collection block 
axial force at stations A, B, C 
segment axial forces at station A 
segment axial forces at station B 
segment axial forces at station C 
axial force function in initialization block 
radial force at stations A, B, C 
X, Y, Z bearing force components 
total axial gap at stations A, B, C 
total axial command at stations A, 

GAAC, GABC, GACC, delayed by DELT 
segment axial gap at station A 
segment axial gap at station B 
segment axial gap at station C 
interval average axial gap at station 
interval minimum axial gap at station 
interval range of axial gap at station A 
interval root mean square axial gap at station A 
interval sum of axial gap at station A 
interval sum of square of axial gap at station A 
interval maximum axial gap at station A 
total radial gap at stations A, B, C 
GRA, GRB, GRC, delayed by DELT 
X, Y, Z rim angular momentum components 
time derivatives of HI, H2, H3 
H1DT, H2DT, H3DT, delayed by DELT 
axial bearing control current at stations A, B, C 
equivalent permanent magnet bias current 
transverse moments of inertia of rim 
polar moment of inertia of rim 

radial bearing control currents at stations A, B, C 
axial controller position and rate gain 
radial controller position and rate gain 


B, C 


A 

A 
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KRB 

KRM 

K1 

M 

NILIM, PILIM 
NP 

0MG1, OMG 2 , OMG 3 
OMG1P, OMG2P, OMG3P 
PI, P2, P3 
P1DT, P2DT, P3DT 
P1DTP, P2DTP, P3DTP 
RAD 
RATIO 
RB 

RBBAR 

RBL 

RBSM 

RBU 

RM 

RPM 

RPMI 

RPMF 

SPA, SPB, SPC 

TEMPI, TEMP 2 

TH1 , TH2 , TH3 

TIME 

TIMEPR 

Tl, T2, T3 

WAA, WAB, WAC 

WAA1 , WAA 2, WAA3, WAA4 

WAB1 , WAB2 , WAB3, WAB 4 

WAC1 , WAC 2, WAC 3, WAC 4 

WAL 

WEIGHT 

WR 

WRL 

X 

XBARO 

XI, X2, X3 

X30 

X30P 

X1DT, X1DT, X3DT 
XI DTP , X2DTP, X3DTP 
Y, YO, YOP 
Z 


radial electromagnet gain 

radial equivalent permanent magnet stiffness 
axial electromagnet gain 
mass of rim 

current limits of power driver 

interval iteration counter 

X, Y, Z rim body rate components 

OMG1 , OMG 2, OMG 3 , delayed by DELT 

X, Y, Z rim linear momentum components 

time derivative of PI, P2, P3 

P1DT, P2DT, P3DT, delayed by DELT 

degree to radian conversion factor 

peak-to-peak axial gap to warp amplitude ratio 

axial displacement of rim center of mass 

average of RB 

interval minimum of RB 

interval sum of RB 

interval maximum of RB 

radius of rim 

rim rotation rate 

initial rpm for simulation 

final rpm for simulation 

position set points for stations A, B, C 

temporary variables in magnetic bearing block 

angular displacement of rim about X, Y, Z axes 

simulation elapsed time 

print interval 

bearing torque components about X, Y, Z axes 

axial warp at stations A, B, C 

segment axial warp at station A 

segment axial warp at station B 

segment axial warp at station C 

axial rim warp amplitude 

rim weight 

radial rim warp at all stations 
radial rim warp amplitude 

temporary variable in initialization block 
equivalent nominal magnetic bearing gap 
X, Y, Z rim displacement components 
initial rigid body axial displacement 
temporary variable in initialization block 
time derivative of XI, X2, X3 
X1DT, X2DT, X3DT, delayed by DELT 
temporary variables in initialization block 
function argument used for FCNY 


AMCD DESCRIPTION 


An AMCD is a spacecraft attitude control effector that provides control 
torques by exchanging spacecraft angular momentum with an internal momentum 
storage device. The storage device is a large spinning rim and the exchange 
mechanism is a set of three noncontacting magnetic bearings which exert forces 
on the spinning rim. A brief description of each of the system components as 
well as the physical layout of the laboratory model follows. 

Orientation of components in the AMCD is shown in figure 1. Three magnetic 
bearing stations, labeled A, B, and C, in the figure, are spaced equidistant 
about the circumference of the rim. At each bearing station, four groups of 
electromagnets are placed about the rim at -14.5, -6, +6, and +14.5 degrees of 
arc relative to the station center. Each group contains two axial and one radial 
bearing element. At each station center, proximity sensors are installed to 
measure the displacement of the rim within the bearing gap. Proximity sensor 
signals and external command signals are input to the controller electronics. 

The controller maintains the rim position by regulating the flow of current in 
the electromagnet coil windings, thereby causing forces to be produced by the 
electromagnets . 


SIMULATION DESCRIPTION 

Figure 2 shows a block diagram of the simulation. The blocks divide the 
simulation into segments according to function. In the following text, the 
function of each of the blocks is described. Appendix A contains a program 
listing to which annotation has been added to enable one to identify the 
sections of code which correspond to the functional blocks. 


Rigid Body and Rim Warp Models 


Models for the rigid body displacements and for axial and radial rim warp 
are contained in the first block. Rim warp is defined as a deformation of the 
rim from an ideal shape due to manufacturing defects or plastic creep due to 
storage and handling. 

Laboratory measurements have indicated that the axial warp may be 
approximated by a sinusoid with two periods per rim revolution (this may be 
visualized as the shape that a large flexible ring would assume if the weight 
were supported by two diametrically opposed points). The amplitude of the 
sinusoid is adjusted to make the model agree with measured data. Axial rim 
warp at each electromagnet and sensor location, 15 points in all, is computed 
using station A center point as zero reference. 

Similar laboratory measurements of the radial rim warp have indicated that 
it may be approximated by a sinusoid with three periods per rim revolution. 

As in the axial case, the amplitude of the sinusoid is adjusted to make the 
model agree with measured data. Radial rim warp is computed at the center of 
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each bearing station using station A as zero reference. Since the bearing 
stations are 120 degrees apart, it may be shown that all the radial rim warp 
components are in phase and, thus, they may be generated from a single model. 

Axial and radial rigid body angular and translational displacements at 
each of the same 15 points mentioned above are derived from figure 1 under the 
assumption of small displacements of the rim in the radial and axial directions. 

Finally, the warp and rigid body contributions to the rim displacements 
are summed for each of the axial and radial components to generate the total 
rim displacement in the bearing gaps. It is this displacement which is 
measured by the proximity sensors and transmitted to the controller. 


Controller Model 

In this simulation, perfect sensors are assumed, therefore, the computed 
bearing gaps are passed directly to the controller and considered to be the 
measured gaps. 

The controller is a proportional plus derivative type which is expressed 
in discrete equation form. Gains are derived from reference 2 and converted to 
equivalent gains for the discrete controller. Simulated proximity sensor 
measurements are combined with set point commands for input to the controller 
equations. Outputs from the controller are voltages which are proportional to 
the required control currents. 

A power amplifier converts the voltage control commands into currents 
which energize each electromagnet circuit. The power amplifiers are modelled 
as having ideal transfer characteristics, with a limit on the maximum available 
output current for that particular bearing segment. The power supply effects 
are modelled as a limit on the total current demanded from all of the power 
amplifiers. After the power amplifier model has processed the command current, 
the total current commanded is computed as the sum of the unsigned magnitudes 
(the sign determines only the direction of the current flow in the electro- 
magnet coil windings). If the upper limit has been reached or exceeded, all 
commands are ratioed to the maximum available current. 


Magnetic Bearing Models 


Each of the twelve nonlinear axial magnetic actuators is an implementation 
of equation 42 in reference 5. An axial magnetic actuator consists of a pair 
of opposed electromagnets, each of which can exert an attractive force. This 
force is proportional to the square of the current flowing through the coil 
windings and inversely proportional to the square of the distance from the rim 
to the magnet pole face. After the forces produced by each of the twelve 
actuators are computed, they are resolved into the rim coordinate system to 
yield two torque components and one force component which will be applied to 
the rigid body dynamics equations. 
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Each of the three radial stations is modelled as a linear system which is 
an implementation of equation 43 in reference 5. A radial station is represented 
by a single electromagnet which exerts an attractive force against the inner 
surface of the rim. Each radial station is opposed by a component of the force 
generated by the other two stations. The forces produced by the three radial 
actuators are resolved into the rim coordinate system to yield two force 
components which will be applied to the rigid body dynamics equations. 


Rigid Body Dynamics 


The rigid body dynamics of the rim are computed by the momentum method as 
described in reference 6. The first step in the computation is to determine 
the rim body rates by dividing the momentum states by the moment of inertia of 
the rim about the respective body axis (the inertia matrix' of the rim is diagonal). 
Then, these body rates are used in the dynamics equation to determine the 
derivative with respect to time of the rim momentum. The momentum derivatives 
are found by solving for the partial derivatives with respect to time from the 
expression for the total derivatives of momentum in the body coordinate system. 
Initially, all cross product terms were incorporated, but subsequent simulation 
data showed that only the terms currently used are significant. 


Data Collection 


In the data collection section of the program, the bearing gap at station 
A, the total required bearing current, and the average displacement of the rim 
along the Z-axis are sampled at each computation interval. Maximum, minimum, 
and average values are computed from these parameters over a time segment defined 
as the print interval. Additionally, the RMS gap, the maximum peak-to-peak 
gap, and a measure of the fraction of the time segment in which the axial 
bearing power amplifier at station A is in a saturated state are computed. 

These data together with the simulation elapsed time and rim RPM are printed at 
the end of each print interval. 


Numerical Integration 


An Adam's second order integrator was selected from a brief comparison 
study of candidate integration schemes wherein speed and error control were 
used as critical selection criteria. At the beginning of the simulation, all 
parameters must be initialized and the steady state conditions of the rim 
dynamics prior to the first computation step must be found. An iterative 
technique is used to solve the rim dynamic and kinematic relationships to yield 
the equilibrium condition. This allows all initial conditions to be established. 
At that point, the simulation begins and proceeds until the final rim RPM value, 
as selected by the user, has been reached. 
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. CONCLUDING REMARKS 


This paper has described a digital computer simulation of a laboratory 
model of an Annular Momentum Control Device (AMCD). The digital computer 
simulation of this device incorporates the following models: six degree-of- 

freedom rigid body dynamics, rim warp, controller dynamics, nonlinear 
distributed element axial bearings, as well as power driver and power supply 
current limits. Annotated FORTRAN IV source code of the computer program is 
included in appendix A. The simulation can be used in the analysis of advanced 
control systems approaches for the laboratory model AMCD. 
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APPENDIX A 


C****************************************************** 

C**** INITIALIZATION BLOCK 1000 ***** 

C****************************************************** 

PROGRAM AMCD 

REAL IAA» IABf IAC» IBIAS» INA1 » INA2* IN A3 
REAL IRA> IRB» IRC 
REAL K1 >KAA»KAR>KRA»KRB»KRMf KRR 
REAL M » NILIM 

DOUBLE PRECISION X * X30» X30P » Y t YO r YOP » Z > FCNY 
LOGICAL FLAG 
1000 PI=4 ♦ *AT AN ( 1 ♦ ) 

RAD=PI/180 * 

DEG=180 ♦ /PI 

C**** SYSTEM PARAMETERS 

M=1 ♦ 54 

WEIGHT ~M*32 . 2 
INA1=5 . 306 
INA2=INA1 
INA3=2 ♦ *INA1 
RM=2 ♦ 625 
KAA-5 ♦ 965E+3 
KAR=2 ♦ 02E+1 
KRA=351 .*12. 

KRR=1 .3*12. 

KRB=1 .55 
KRM=170 . *12 ♦ 

Kl=3 ♦ 4476E-3/144 . 0/16 . 0 
PILIM=13. 33 
NILIM=-PILIM 
CRLIM-30 . 

IBIAS=14. 96999297 
XBARO=. 14/12. 

SPA=. 025/12. 

SPB=SPA 

SPC=SPA 

WRL=.002/12. 

C**** INPUT CASE DATA 

1100 WRITE ( 5 1 1200 ) 

1200 FORMAT ( ' SELECT DELTA T (END PROGRAM IF NEGATIVE)') 

ACCEPT 1300 » BELT 
1300 FORMAT ( F12 ♦ 3 ) 

IF ( BELT ♦ LE ♦ 0 ♦ ) GOTO 9999 
WRITE ( 5 r 1400 ) 

1400 FORMAT ( ' SELECT CASE PARAMETERS ...V 

1 # PRINT INTERVAL 1 1NIT SPIN RPM» FINAL SPIN RPM»WARP AMPLITUDE') 
ACCEPT 1300 > DELTPR > RPMI » RPMF » WAL 
C**** DETERMINE STEADY STATE BEARING GAP 

FCNY(Z)~12.*K1*((( IBIAS+KAA* < SPA-Z ) ) / ( XBARO-Z ) )**2 
1 -< (IBIAS-KAA*(SPA-Z))/(XBARO+Z) ) **2 ) -WEIGHT 
X30P= ♦ 010/12 ♦ 

YOP=FCNY ( X30P ) 

X30=. 020/12. 

Y0=FCNY(X30) 

1500 X= ( X30P*Y0-X30*Y0P ) /( YO-YOP ) 

Y=FCNY ( X ) 



c**** 


X30P=X30 

YOF-YO 

X30=X 

YO=Y 

IF ( ABS ( Y ) .GT ♦ 1 ♦ E-10 ) GOTO 1500 
TH3=0 ♦ 

WAA=WAL*SIN(2* *TH3 ) 

WAB=WAL*<SIN(2.*TH3>*C0S<240.*RAD)+C0S(24*TH3>*SIN<240.*RAn) > 

WAC=WAL*<SIN<2,*TH3>*C0S<120.*RAri)+C0S<2.*TH3>*SIN(120.*RAr.i) ) 

GAA=-UAA+X30 

GAB=-WAB+X30 

GAC=-WAC+X30 

TH1=(GAA-GAB)/(RM*SQRT(3» ) ) 

TH2=<2.*GAC-<GAA+GAB> >/<RM*3. > 

X1=0. 

X2=0 . 

X3= ( GAA+GAB+GAC ) /3 ♦ 

GAACF-X30-SPA 

GABCP=X30-SF'B 

GACCP=X30-SFC 

UR=WRL*SIN<3.*TH3> 

GRA=X1*SIN<30.*RAIU+X2*C0S<30.*RAD)+UR 

GRB=Xl*SIN<30.*RAD)-X2*C0S(30.*RAn>+WR 

GRC=-X1+WR 

GRAP=GRA 

GRBF'=GRB 

GRCP-GRC 

COMPUTE INITIAL CONDITIONS 
TIME=0 ♦ 

TIMEPR--DELT/2* 

OMGl^O ♦ 

0MG2=0 ♦ 

0MG3=RPMI*2.*PI/60. 

X 1 DT =0 « 

X2DT =0 ♦ 

X3DT -0 ♦ 

H1=0MG1*INA1 

H2-0MG2*INA2 

H3=0MG3*INA3 

P1«X1DT*M 

F'2=X2DT*M 

P3=X3DT*M 

T 1 =0 ♦ 

T2=0 ♦ 

T3=l.l 

H1DT=T1-0MG2*H3 

H2DT=T2+0MG1*H3 

H3DT=T3-0MG1*H2+0MG2*H1 

P1DT=0. 

P2DT =0 . 

P3DT -0 ♦ 

H1DTP-H1DT 

H2DTP-H2DT 

H3DTP=H3DT 

P1DTF-P1DT 

P2DTP=P2DT 

F'3DTF’=P3DT 

0MG1F'~0MG1 

0MG2P-0MG2 
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c**** 


0MG3F'=QMG3 

X1DTP=X1DT 

X2DTP=X2DT 

X3DTP=X3DT 

CLEAR STATISTICAL VARIABLES 
GPSM=0 . 

GPSQSM=0 ♦ 

GPL=X30 
GPU=GPL 
CRSM=0 ♦ 

CRL=KAA*(ABS(GAACP>+ABS<GABCP)+ABS(GACCP> ) 

CRU=CRL 

RBSM=0 ♦ 

RBL= ( GAACP+GABCP+GACCP ) /3 ♦ 

RBU=RBL 
CNG=0 . 

CPS=0. 

NP=0 

FLAG= ♦ FALSE . 


C****************************************************** 

C**** RIGID BODY AND RIM WARP BLOCK 2000 ***** 

C************************************ ****************** 


c**** 

2000 


C**** 


AXIAL WARP 
CONTINUE 

WAA1=WAL*(SIN(2« *TH3) *COS< -29 «*RAD)+C0S(2»*TH3)*SIN( 
WAA2=WAL*<SIN(2* *TH3 ) *COS ( -12 ♦ *RAD)+CGS(2**TH3)*SIN< 
WAA = WAL* SIN ( 2 ♦ *TH3) 

WAA3=WAL*<SIN(2. *TH3> *COS( 12 . *RAD ) + COS < 2 . *TH3 ) *SIN ( 
WAA4=WAL*(SIN(2«*TH3) *COS ( 29 ♦ *RAD ) +COS ( 2 . *TH3 ) *SIN< 
WAB1=WAL*(SIN(2* *TH3 ) *COS <211«*RAD)+C0S(2**TH3)*SIN( 
WAB2=WAL* ( SI N ( 2 . *TH3 ) *COS ( 22B . *RAD ) +COS ( 2 . *TH3 ) *SIN ( 
WAB =WAL*(SIN(2* *TH3 )*COS (240«*RAD)+C0S(2**TH3)*SIN( 
WAB3=WAL*(SIN(2« *TH3 )*COS( 252 . *RAD )+C0S(2»*TH3)*SIN( 
WAB4=WAL*(SIN(2* *TH3 ) *COS ( 269 ♦ *RAD )+C0S(2«*TH3)*SIN( 
WAC1=WAL*(SIN(2.*TH3) *COS ( 91 ♦ *RAD ) +COS ( 2 ♦ *TH3 ) *SIN ( 
WAC2=WAL*(SIN(2« *TH3 )*COS( 108 . *RAD ) +COS ( 2 ♦ *TH3 ) *SIN ( 
WAC =WAL*(SIN(2. *TH3 > *COS ( 120 . *RAD ) +COS (2«*TH3)*SIN( 
WAC3=WAL*(SIN(2« *TH3 ) *COS( 132 « *RAD ) +COS ( 2 . *TH3 ) *SIN ( 
WAC4=WAL*(SIN(2«*TH3) *COS <149«*RAD)+C0S(2»*TH3)*SIN( 
RADIAL WARP 
WR~WRL*SIN ( 3 . *TH3 > 


-29 ♦ *RAD ) ) 
-12 ♦ *RAD ) ) 

12 ♦ *RAD ) ) 
29 . *RAD > ) 
211 ♦ *RAD ) ) 
228 ♦ *RAD ) ) 
240 ♦ *RAD ) ) 
252 ♦ *RAD ) ) 
269 « *RAD ) ) 
91 . *RAD ) ) 
108 ♦ *RAD ) ) 
120 . *RAD ) ) 
132 ♦ *RAD ) ) 
1 49 » *RAD ) ) 


C**** AXIAL BEARING GAPS 

GAA1=RM*<TH1*C0S< 15 ♦ 5*RAD > -TH2*SIN ( 15.5*RAD> HX3+WAA1 
GAA2-RM* ( TM1 *COS ( 24 . *RAD ) -TH2*SIN ( 24 . *RAD ) ) +X3+WAA2 
GAA=RM*(TH1*C0S(30» *RAD > -TH2*SI N ( 30 « *RAD ) )+X3+WAA 
GAA3=RM* ( TH1*C0S ( 36 . *RAD ) -TH2*SIN < 36 ♦ *RAD ) ) +X3+WAA3 
GAA4=RM* ( TH1 *COS ( 44 . 5*RAD ) -TH2*SIN ( 44 , 5*RAD ) ) +X3+WAA4 
GAB1=RM*(-TH1*C0S(44»5*RAD)-TH2*SIN(44.5*RAD) M-X3+WAB1 
GAB2-RM* ( -TH1 *COS ( 36 . *RAD ) -TH2*SIN ( 36 . *RAD ) ) +X3+WAB2 
GAB=RM*(“TH1*C0S(30. *RAD ) -TH2*SIN (30«*RAD) )+X3+WAB 
GAB3=RM* ( -TH1 *COS ( 24 . *RAD ) -TH2*SIN ( 24 . *RAD ) ) +X3+WAB3 
GAB4=RM* ( -TH1 *COS ( 1 5 ♦ 5*RAD ) -TH2*SIN ( 15 . 5*RAD ) ) +X3+WAB4 
GAC1=RM*(~TH1*SIN(14* 5*RAD ) +TH2*C0S ( 14 * 5*RAD ) H-X3+WAC1 
GAC2=RM*(“TH1*SIN(6* *RAD ) +TH2*C0S < 6 ♦ *RAD) )+X3 + WAC2 
GAC=RM*TH2+X3+WAC 

GAC3=RM* ( TH1 *SI N ( 6 ♦ *RAD ) +TH2*C0S ( 6 ♦ *RAD ) ) +X3+WAC3 
GAC4=RM*(TH1*SIN(14« 5*RAD ) +TH2*C0S ( 14 ♦ 5*RAD ) )+X3 + WAC4 
C**** RADIAL BEARING GAPS 

GRA=X1*SIN(30.*RAD) +X2*C0S ( 30 ♦ *RAD ) +WR 
GRB=X 1 *SI N ( 30 ♦ *RAD ) -X2*C0S ( 30 ♦ *RAD ) +WR 
GRC=-X1+WR 
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0 * ***************************************************** 
C**** BEARING CONTROLLER BLOCK 3000 ***** 

i“' 'Aj 'X* X* X- X» *,L X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 

^ ^ *• •* ■ ** «p ^ •* ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^p ^p -P p p p p p p p p p p p p p p p p p p p p p 

C**** AXIAL CONTROLLER 

GAAC~6AA-~SPA 

I AA = ~KAA*GAAC--KAR* (GAAC-GAACP ) /BELT 

GAACP=GAAC 

GABC-GAB-SPB 

I AB=-K A A*G ABC-K AR* ( GABC-GABCP ) /BELT 

GABCP=GABC 

GACC-GAC--SPC 

I AC=~KAA*GACC~KAR* ( GACC-GACCP ) /BELT 
GACCP^GACC 

C**** RADIAL CONTROLLER 

IRA=-KRA*GRA-KRR*(GRA-GRAP)/DELT 


C**** 

3100 


3200 


3300 

3400 

3500 


3600 

3700 

err** 


3800 


GRAP-GRA 

I RB~~KRA*GRB-KRR* ( GRB-GRBP ) /BELT 
GRBP=GRB 

IRC=~KRA*GRC-KRR* ( GRC-GRCP) /BELT 
GRCP=GRC 

CURRENT LIMITER 


I F ( I A A . GT . N I L I M ) 
CNG-CNG+1 . 

I AA-NIL.IM 
GOTO 3300 

GOTO 

3200 

IF < I AA ♦ LT . PILIM ) 
CPS=CPS+1 . 
IAA“PILIM 

GOTO 

3300 

I F < I AB ♦ GT ♦ NI L I M ) 
TAB--NILIM 
GOTO 3500 

GOTO 

3400 

IF < I AB * L.T .PILIM ) 
I AB=PILIM 

GOTO 

3500 

I F ( I AC * GT . N I L I M ) 
I AC-NILIM 
GOTO 3700 

GOTO 

3600 

TF( I.AC.LT* PILIM) 

IAC-PILIM 

CONTINUE 

GOTO 

3700 


POWER SUPPLY LIMITER 

»: : R = : ABS ( I A A > TABS ( I AB)+ABS ( I AC) TABS ( IRA ) TABS ( IRB ) TABS (IRC) 

IF'CR.LT.CRLIM) GOTO 3800 

I A A = I A A * C R L I M / C R 

IAB"IAB#CRLIM/CR 

I A C - 1 A C * C R L I M / C R 

IRA~IRA*CRL.IM/CR 

IRB=IRB#CRLIM/CR 

I R C ™ I R C * C R L. 1 M / C R 

CR=CRLIM 

CONTINUE 
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*:■ * ***** W********************************************** 

»'•'**** MAGNETIC BEARING MODEL BLOCK 4000 ***** 

v +***************************************************** 

TEMP 1=( IBIAS+IAA>/( XBARO-GAA 1 ) 

TEMF‘2= < I B I AS- I AA ) / ( XBAR0+GAA1 ) 

FA 1 -K 1 * < TEMP 1 *TEMP1 -TEMP2*TEMP2 ) 

TEMPI ~ ( IBI AS+IAA) / ( XBAR0-GAA2 ) 
f EMP2- ( IBI AS-I AA ) / ( XBAR0+GAA2 ) 
FA2=K1*<TEMP1*TEMP1-TEMP2*TEMP2> 

TEMPI - ( IBI AS + IAA ) / ( XBAR0-GAA3 ) 

TEMP2* (IBI AS-I AA >/ ( XBAR0+GAA3 > 

FA3-K1* < TEMPI *TEMF’1-TEMP2*TEMP2 ) 

TEMPI =< IBI AS+IAA ) / ( XBAR0-GAA4 ) 

TEMP2- ( IBI AS-I AA > / ( XBAR0+GAA4 > 

FA4-K1* ( TEMP1*TEMP1-TEMP2*TEMP2 ) 

FA-FA 1+FA2+FA3+FA4 

TEMP 1 - ( I B I AS+ 1 AB ) / < XBARO-GAB 1 ) 

TE MP2 - <1 B I A S - 1 AB ) / ( XB ARO f G AB 1 ) 

FB I. =K1*( TEMP 1*TEMP1-TEMP2*TEMP2> 

TEMP I- ( IBIAS+IAB ) / ( XBAR0-GAB2 ) 

TEMP2- < IE<IAS-I AB ) / ( XBAR0+GAB2) 

FB2-K 1 * ( TEMPI *TEMP1-TEMP2*TEMP2 ) 

TEMP 1 = < IBIAS+IAB ) / ( XBAR0-GAB3 ;■ 

TEMP2= < IBI AS-IAB ) / ( XBAR0+GAB3 ) 

F r B3-Kl* ( TEMP1*TEMP1-TEMP2*TEMP2 ) 

TEMPI =( IBIAS+IAB > / < XBAR0-GAB4 ) 

TEMP2= ( IBI AS-I AB) / (XBAR0+GAB4 ) 

FB4-K1* ( TEMP1*TEMP1-TEMP2*TEMP2 ) 

FB-FB1+FB2+FB3+FB4 

TEMP 1 = < I B I AS+ I AC ) / ( XBARO-GAC 1 ) 

TE"MF : '2- ( IBI AS-I AC) / ( XBAR0+GAC1 ) 

F C 1 - K 1 * < T E M P 1 * T E M P 1 - T E M P 2 * T E M P 2 ) 

TEMP 1 » ( I B I AST I AC ) / < XBAR0-GAC2 ) 

TEMP2- ( 1 B I AS- I AC ) / ( XBAR0+GAC2 ) 

F C 2 * K 1 * <T E M P 1 * T E M P 1 - T E M P 2 * T E M p 2 ) 

TEMP 1 - < I BIAB+ I AC ) / ( XBAR0-GAC3 > 

TEMP2~( I BIAS- I AC)/ (XBAR0+GAC3) 

FC3=K 1 * < TEMP 1 *TEMP1-TEMP2*TEMP2 ) 

TEMP 1 ~ ( I B I AS + I AC ) / ( XBAR0-GAC4 ) 

TEMP 2 = (IBI AS-I AC ) / ( XBAR0+GAC4 ) 

F ; 'C4~K1 * ' TEMPI *TEMP1 -TEMP2*TEMP2 ) 

FC=FClfFC2+FC3+FC4 

**+* TRANSFORM BEARING FORCES INTO RIM SYSTEM 

T l -'EM* < FA 1 -FB4 ) *COS ( 15 * 5*RAB ) < FA2-FB3 ) *C03 < 24 . *RA r J ) 
.1. i ( FA3-FB2 ) *COS ( 36 . *RAD ) + ( FA4-F r Bl ) *COS < 44 . 5*RAD > 

- > • F C 3 - F ' C 2 ) * S I N ( 6 . * R A D ) + ( F C 4 - F C 1 ) * S I N <1 4 . 5 * R A D ) > 

T2 ~RM* ( ■ FC2+FC3 ) *COS ( 6 . *RAD ) + ( FC 1 f FC4 > *COS (14* 5*RAD > 

1 - ( FA4+FB1 ) * S I N ( 44 * 5*RAD) - ( FA3+FB2 ) *SIN ( 36 s *RAIO 

2 - ( FA2+F B3 ) *S I N ( 24 . *RAD ) - ( F A 1 +FB4 ) *S I N ( .1 5 ♦ 5* RAD ) ) 
F3-FA+FB+FC 

IF** RADIAL BEARINGS 

FRA-KRB* I RA+ KRM*GRA 
F R B = K R B * I R B + K R M * G R B 
FROKRB* I RC+KRM*GRC 
F ; 1 = ( F R A t F F i.i ) * S I N < 3 0 . * R A D > - F R C 
F2 ::: < FRA FRB > *COS ( 30 . *RAB ) 
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I*' ^ ^ ^ ^ ^ ^ ^ Xf <j» Jj ^ ^ fti« ^ lL ^ ^ ^ ||j ^ X* «L t ja J* at* ala til «|j yK tli ^1* ala j l. ala tlf » l| ill ala Ji «!■ ala i|a a^ 

w ^ ^ ^ip t' ^ ^ * T^ * * ^ ^ ^ ^ * * * * * ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

C**** RIM DYNAMICS BLOCK 5000 ***** 

r' ^ \l^ «Xr »X ^ ^ ^ Xf ^ Xr ^ Xf Xf Xf ^ llf> Xf ^ ^ ^ ^ ^ ^ ^ ^ Uj ^ ^ J« ll » ^Ji Ji all it. J. afa ala J. J. ala >K aja ill ^a ala at* «i« ala ala 

W ^ ^ ^P ^ ^p ^ MP * ^ ^ ^ ^ ^ ^p * ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

0MG1=H1/INA1 

0MG2~H2/INA2 

0MG3=H3/INA3 

X1DT=P1/M 

X2DT=P2/M 

X3DT=P3/M 

H1DT=T1-0MG2*H3 

H2DT=T2+0MG1*H3 

H3DT=T3-0MG1*H2+0MG2*H1 

P1DT =F 1 

P2DT =F2 

P3DT =F3 -WEIGHT 


c ****************************************************** 

C**** DATA COLLECTION AND I/O BLOCK 6000 ***** 

I ^Xr ^ ^Lr ^ ^ ^ ^ \Xf ^ ^ iX» ^ ^ ^ JLp ^ ^ ^ ^ ^ ^ ^ ^ ^ *U «l» ^ »>• *|j ^ «!■ «l» ^ «Ij tL tf* ^ U* y* J. »!> •!* ^i* 

w ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

GF‘SM=GPSM+GAA 
GPSGSM=GF'SGSM+GAA*GAA 
IF < GAA ♦ GT . GPU ) GPU=GAA 
IF < GAA ♦ LT « GPL ) GPL=GAA 
CRSM=CRSM+CR 
IF < CR ♦ GT • CRU ) CRU=CR 
IF ( CR « LT ♦ CRL ) CRL=CR 
RB= ( GAA+GAB+GAC ) /3 ♦ 

RBSM=RBSM+RB 
IF ( RB . GT ♦ RBU ) RBU=RB 
IF (RBtLT .RBL) RBL=RB 
NF'=NF'+ 1 

C**** PREPARE DATA FOR OUTPUT 

IF ( TIME . LT ♦ TIMEF'R > GOTO 6300 
TIMEPR=TIMEPR+DELTPR 
RF‘M=30 « *0MG3/F'I 



WRITE ( 6 7 6000 ) 

TIME»RPM 




6000 

FORMAT ( 7 TIME 

='F10.3»2Xf 7 RPM 

7 F10 ♦ 3 ) 




IF (FLAG) GOTO 6100 





WRITE ( 6 7 6010 ) 

H1DT »H2DT »H3DT »P3DT 




6010 

FORMAT ( 7 HI DT 

=' 1PE10.3»2X, 7 H2DT 

= 7 El 0 « 3 r 2X t 7 H3DT 

= 7 E10.3f 

2X7 


1 7 P3DT 

= ' E10 . 3 ) 





WRITE ( 6 f 6020 ) 

HI »H2>H3»P3 




6020 

FORMAT ( ' HI 

= 7 1PE10 ♦ 3 » 2X r 7 H2 

= 7 E10 ♦ 3 » 2X f 7 H3 

= 7 E 1 0 ♦ 3 7 

2X7 


1 7 F'3 

=/ E10«3> 





WRITE ( 6 1 6030 ) 

0MG1 »0MG2»0MG3»X3DT 




6030 

FORMAT ( 7 0MG1 

= ' 1F*E10.3»2X» 7 0MG2 

= 7 E10*3»2X» 7 0MG3 

= 7 E10 4 3 7 

2X7 


1 7 X3DT 

= 7 E10.3> 





WRITE ( 6 t 6040 ) 

TH1 1 TH2» TH3 7 X3 




60-40 

FORMAT < 7 TH1 

= 7 1PE10 ♦ 3 1 2X > 7 TH2 

= 7 E 1 0 ♦ 3 » 2X r 7 TH3 

= 7 E 1 0 4 3 7 

2 X 7 


1 7 X3 

= 7 E10.3> 





WRITE ( 6 r 6050 ) 

T 1 7 T2 7 F3 




6050 

FORMAT <' T 1 

= 7 lPE10.3r2Xr 7 T2 

- 'E10.3f 2X r 7 F3 

= 7 E 1 0 4 3 ) 



WRITE ( 6 7 6060 ) 

FA 7 FB 7 FC 




6060 

FORMAT (' FA 

= 7 1F'E10.3»2X» 7 FB 

= 7 E10 ♦ 3 ? 2X > 7 FC 

= 7 E 1 0 4 3 ) 



WRITE (6t 6070 ) 

IAA 7 IAB 7 IAC 




6070 

FORMAT ( 7 IAA 

= 7 1PE10.3j 2X r 7 1 AB 

= 7 E 1 0 . 3 » 2X » 7 1 AC 

= 7 E10 4 3 ) 


C**** 

COMPUTE GAP STATISTICS 




6100 

DI 01=0 ♦ 






IF ( NP ♦ GT < 0 ) DIV1«1./NP 
GF‘BAR=GF’SM*DIV1 
CRBAR=CRSM*DI VI 
RBBAR=RBSM*DIV1 
DI V2-0 « 

IF(NP.GT.l) DIV2=1./(NP-1> 

GPRMS=SQRT < ( GPSGSM~NP*GPBAR*GPBAR ) *DI 02 ) 

GPBAR=12 . *GF'BAR 

GPRMS=12.*GPRMS 

GPU=12.*GPU 

GPL=12.*GPL 

RBU= 1 2 . *RBU 

RBL=12 . *RBL 

RBBAR=12**RBBAR 

GPPF-ABS(GPU--GPL) 

CNG=100 « *CNG*DIV1 
CF'S=100 • *CF‘S*DIV1 
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RATI0=GPPP/<WAL*24. ) 

WRITE < 6 > 6200 ) GPBAR > GPRHS * GPL » GPU 
6200 FORMAT ( 7 GPBAR = 7 1PE10*3> 2X» 7 GPRMS = ' E10 ♦ 3 r 2X r 7 GPL 
1 7 GPU ="E10*3) 

WR I TE ( 6 » 621 0 > GPPP » CRBAR , CRL » CRU 
6210 FORMAT < ' GPPP = 7 1PE10«3»2X* 'CRBAR ='E10.3>2X» 'CRl. 
1 7 CRU = 7 E10«3) 

WRITE < 6 > 6220 > CNG r CPS f RBL » RBU 
6220 FORMAT ( 7 CNG = 7 1PF10«3»2X* 'CPS = 7 F10 ♦ 3 > 2X , 7 RBL 
1 7 RBU = 7 E10»3) 

GPSM=0 . 

GPSQSM=0. 

GPL=GPBAR/12 ♦ 

GPU=GPL 
CRSM=0 « 

CRL=CRBAR 
CRU=CRL 
RBSM=0 ♦ 

RBL=RBBAR/12 . 

RBU=RBL 
CNG=0 ♦ 

CPS=0. 

NP~0 

FLAG= , TRUE . 

6300 CONTINUE 


= 7 E10 . 3 > 2X 


= 7 E10 ♦ 3 > 2X 


= 7 E10 ♦ 3 > 2X 



C****************************************************** 
C**** NUMERICAL INTEGRATION BLOCK 7000 ***** 

^1/ ^ i±r ^ ^ ^ ^ \Lr ^ tb >b sb \b ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ tb ^b ^b >l» d* ^b tb <b J# ^b b t b ^b «b d* ib ^b d« J. ^b J. 

»/ * * nN * ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

H1=H1 + DELT* ( 3 « *HlDT-HlDTF')/2 ♦ 
H2=H2+DELT*(3.*H2DT-H2DTP)/2. 

H3=H3+DELT* < 3 . *H3DT-H3DTP ) /2 . 
Pl=Pl+DELT*<3.*PlDT-PlDTP)/2. 

P2=P2+DELT* ( 3 . *P2DT-P2DTP ) /2 ♦ 

P3-P3+DELT* ( 3 . *P3DT-P3DTP ) /2 ♦ 

TH1=TH1+DELT*(3. *0MG1-0MG1P) /2 ♦ 

TH2=TH2+DELT* ( 3 # *0MG2-0MG2P ) /2 ♦ 
TH3=TH3+DELT*(3.*OMG3-OMG3P>/2. 

X1=X1+DELT* ( 3 ♦ *XlDT-XlDTP>/2 . 

X2=X2+DELT* ( 3 ♦ *X2DT-X2DTP)/2 . 

X3=X3+DELT* ( 3 . *X3DT-X3DTP ) /2 . 

H1DTP=H1DT 

H2DTP=H2DT 

H3DTP=H3DT 

P1DTP=P1DT 

P2DTP=P2DT 

P3DTP=P3DT 

0MG1P=0MG1 

0MG2P=0MG2 

0MG3P=0MG3 

X1DTP=X1DT 

X2DTP=X2DT 

X3DTP=X3DT 

IF ( TH3 . GT . 360 ♦ *RAD ) TH3=TH3-360 . *RAD 

TIME=TIME+DELT 

IF ( RF'M . LT ♦ RPMF ) GOTO 2000 

GOTO 1100 

C*** TERMINATE PROGRAM 

9999 CONTINUE 
END 


18 














1. Report No. 2. Government Accession No. 

NASA TM-81797 

3. Recipient's Catalog No. 

4. Title and Subtitle DESCRIPTION OF A 

DIGITAL COMPUTER SIMULATION OF AN ANNULAR MOMENTUM 

5. Report Date 
February 1981 

CONTROL DEVICE (AMCD) LABORATORY TEST MODEL 

6. Performing Organization Code 

7. Author(s) 

8. Performing Organization Report No. 

Charles T. Woolley and Nelson J. Groom 



10. Work Unit No. 

9. Performing Organization Name and Address 

506-54-83-01-00 

NASA Langley Research Center 
Hampton, Virginia 23665 

11. Contract or Grant No. 


13. Type of Report and Period Covered 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Waslrrngton, DC 20546 

Technical Memorandum 

14. Army Project No. 

15. Supplementary Notes 


6. Abstract 


A description of a digital computer simulation of an Annular Momentum Control 
Device (AMCD) laboratory model is presented. The AMCD is a momentum exchange device 
which is under development as an advanced control effector for spacecraft attitude 
control systems. 

The digital computer simulation of this device incorporates the following 
models: six degree-of-freedom rigid body dynamics, rim warp, controller dynamics, 

non-linear distributed element axial bearings, as well as power driver and power 
supply current limits. An annotated FORTRAN IV source code listing of the computer 
program is included in appendix A. 


17. Key Words {Suggested by Authors) ) 18. Distribution Statement 

Digital computer simulation 

Annular momentum control device Unclassified - Unlimited 

Magnetic bearings 

Subject Category 61 


19. Security Classif. (of this report} 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price* 

Unclassified , 

Unclassified 

20 

A0 2 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


















